module phys_debug

! This module contains subroutines that are intended to print diagnostic
! information about a single column to the log file.  The module uses the
! phys_debug_util module to locate the column using local indices of the
! physics decomposition.  The reason to encapsulate the debug print
! statements into a subroutine is to minimize the impact of the debugging
! code on the development code.  This makes it easier to maintain debug
! code on a branch which can be updated to the head of the trunk when
! convenient.  Minimizing the footprint of debug code helps to minimize the
! chances of conflicts that can occur while updating between the debug code
! and changes that have occurred on the trunk.

! This module also contains some generic code to write a collection of variables
! from the physics state object to the history file.


use shr_kind_mod,    only: r8 => shr_kind_r8
use phys_debug_util, only: phys_debug_col
use ppgrid,          only: pcols, pver, pverp
use physics_types,   only: physics_state, physics_ptend
use camsrfexch,      only: cam_in_t
use cam_logfile,     only: iulog
use cam_history,     only: addfld, add_default, outfld
use constituents,    only: pcnst, cnst_get_ind, cnst_name
use rad_constituents,only: rad_cnst_get_info, rad_cnst_get_mode_num_idx, rad_cnst_get_mam_mmr_idx

implicit none
private
save

! Generic routines for writing physics state to history
public :: &
   phys_debug_state_init, &
   phys_debug_state_out


! The following routines are intended as examples for printing various diagnostic
! information from a single column
public :: &
   phys_debug_vdiff1,   &
   phys_debug_shallow1, &
   phys_debug_strat1,   &
   phys_debug_srf1,     &
   phys_debug_hbdiff1,  &
   phys_debug_flux1,    &
   phys_debug_flux2

! constituent indices
integer :: ixcldliq, ixcldice, ixnumliq, ixnumice
! logical array to identify constituents that are mode number concentrations
logical :: cnst_is_mam_num(PCNST)
! logical array to identify constituents that are mode specie mass mixing ratios
logical :: cnst_is_mam_mmr(PCNST)


!===================================================================================================
contains
!===================================================================================================

subroutine phys_debug_state_init(tags)

   ! Write selected state variables from locations corresponding to the 'tags'.
   ! This method does the addfld/add_default calls.

   ! ** Note ** This init needs to be called after the rad_constituent module has finished
   !            its init.


   character(len=8), dimension(:), intent(in) :: tags

   ! local variables
   integer :: i, icnst, j, nmodes, nspec, ntags
   character(len=5) :: unit
   !-----------------------------------------------------------------------------------------------

   ! Initialize some constituent indices and names
   call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
   call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
   call cnst_get_ind('NUMLIQ', ixnumliq, abort=.false.)
   call cnst_get_ind('NUMICE', ixnumice, abort=.false.)

   ! Set arrays to identify the modal aerosol constituents
   cnst_is_mam_num = .false.
   cnst_is_mam_mmr = .false.
   call rad_cnst_get_info(0, nmodes=nmodes)
   do i = 1, nmodes
      call rad_cnst_get_mode_num_idx(i, icnst)
      cnst_is_mam_num(icnst) = .true.
      call rad_cnst_get_info(0, i, nspec=nspec)
      do j = 1, nspec
         call rad_cnst_get_mam_mmr_idx(i, j, icnst)
         cnst_is_mam_mmr(icnst) = .true.
      end do
   end do

   ! For each supplied tag, issue addfld/add_default calls for the state variables
   ! with the tag name appended.
   ntags = size(tags)
   do i = 1, ntags

      call addfld('T_'//trim(tags(i)), (/ 'lev' /), 'A', 'K', 'T at tag '//trim(tags(i)))
      call add_default('T_'//trim(tags(i)), 1, ' ')

      call addfld('U_'//trim(tags(i)), (/ 'lev' /), 'A', 'm/s', 'U at tag '//trim(tags(i)))
      call add_default('U_'//trim(tags(i)), 1, ' ')

      call addfld('V_'//trim(tags(i)), (/ 'lev' /), 'A', 'm/s', 'V at tag '//trim(tags(i)))
      call add_default('V_'//trim(tags(i)), 1, ' ')

      call addfld('QV_'//trim(tags(i)), (/ 'lev' /), 'A', 'kg/kg', 'QV at tag '//trim(tags(i)))
      call add_default('QV_'//trim(tags(i)), 1, ' ')

      if (ixcldliq > 0) then
         call addfld('QL_'//trim(tags(i)), (/ 'lev' /), 'A', 'kg/kg', 'QL at tag '//trim(tags(i)))
         call add_default('QL_'//trim(tags(i)), 1, ' ')
      end if

      if (ixcldice > 0) then
         call addfld('QI_'//trim(tags(i)), (/ 'lev' /), 'A', 'kg/kg', 'QI at tag '//trim(tags(i)))
         call add_default('QI_'//trim(tags(i)), 1, ' ')
      end if

      if (ixnumliq > 0) then
         call addfld('NL_'//trim(tags(i)), (/ 'lev' /), 'A', '1/kg', 'NL at tag '//trim(tags(i)))
         call add_default('NL_'//trim(tags(i)), 1, ' ')
      end if

      if (ixnumice > 0) then
         call addfld('NI_'//trim(tags(i)), (/ 'lev' /), 'A', '1/kg', 'NI at tag '//trim(tags(i)))
         call add_default('NI_'//trim(tags(i)), 1, ' ')
      end if

      do j = 1, pcnst

         if (cnst_is_mam_num(j) .or. cnst_is_mam_mmr(j)) then

            unit = '1/kg'
            if (cnst_is_mam_mmr(j)) unit = 'kg/kg'

            call addfld(trim(cnst_name(j))//'_'//trim(tags(i)), (/ 'lev' /), &
               'A', trim(unit), trim(cnst_name(j))//' at tag '//trim(tags(i)))
            call add_default(trim(cnst_name(j))//'_'//trim(tags(i)), 1, ' ')

         end if
      end do

   end do

end subroutine phys_debug_state_init

!===================================================================================================

subroutine phys_debug_state_out(tag, state)

   ! Dump selected state variables at location corresponding to given tag.

   character(len=*),    intent(in) :: tag
   type(physics_state), intent(in) :: state

   ! Local variables
   integer :: lchnk, j
   !-----------------------------------------------------------------------------------------------

   lchnk = state%lchnk

   call outfld('T_'//trim(tag),  state%t,        pcols, lchnk)
   call outfld('U_'//trim(tag),  state%u,        pcols, lchnk)
   call outfld('V_'//trim(tag),  state%v,        pcols, lchnk)
   call outfld('QV_'//trim(tag), state%q(:,:,1), pcols, lchnk)
   if (ixcldliq > 0) call outfld('QL_'//trim(tag), state%q(:,:,ixcldliq), pcols, lchnk)
   if (ixcldice > 0) call outfld('QI_'//trim(tag), state%q(:,:,ixcldice), pcols, lchnk)
   if (ixnumliq > 0) call outfld('NL_'//trim(tag), state%q(:,:,ixnumliq), pcols, lchnk)
   if (ixnumice > 0) call outfld('NI_'//trim(tag), state%q(:,:,ixnumice), pcols, lchnk)

   do j = 1, pcnst
      if (cnst_is_mam_num(j) .or. cnst_is_mam_mmr(j)) then
         call outfld(trim(cnst_name(j))//'_'//trim(tag), state%q(:,:,j), pcols, lchnk)
      end if
   end do

end subroutine phys_debug_state_out


!===================================================================================================

subroutine phys_debug_vdiff1(state, kvm, tautotx, tautoty)

   type(physics_state), intent(in) :: state             ! Physics state variables
   real(r8),            intent(in) :: kvm(pcols,pverp)  ! eddy diffusivity for momentum [m2/s]
   real(r8),            intent(in) :: tautotx(pcols)    ! u component of total surface stress
   real(r8),            intent(in) :: tautoty(pcols)    ! v component of total surface stress

   integer :: icol
   !-----------------------------------------------------------------------------

   icol = phys_debug_col(state%lchnk)
   if (icol > 0) then
      write(iulog,*) ' vert_diff: kvm, u ', kvm(icol,pver), kvm(icol,pver-1), &
         state%u(icol, pver), state%u(icol,pver-1), &
         state%v(icol, pver), state%v(icol,pver-1)
      write(iulog,*) ' vert_diff: taux, tauy ', tautotx(icol), tautoty(icol)
   endif

end subroutine phys_debug_vdiff1

!================================================================================

subroutine phys_debug_shallow1(state, ptend, nstep, prec_cmf, rliq2, ztodt, kmx)

   use constituents, only: cnst_get_ind

   type(physics_state), intent(in) :: state             ! Physics state variables
   type(physics_ptend), intent(in) :: ptend             ! Physics process tendencies
   integer,             intent(in) :: nstep
   real(r8),            intent(in) :: prec_cmf(pcols)   ! total precipitation from Hack convection
   real(r8),            intent(in) :: rliq2(pcols)      ! vertical integral of liquid from shallow scheme
   real(r8),            intent(in) :: ztodt             ! physics time step
   integer,            intent(out) :: kmx

   integer  :: icol, k, kmn, ixcldliq
   real(r8) :: qtmx, qtmn
   !-----------------------------------------------------------------------------

   icol = phys_debug_col(state%lchnk)
   if (icol > 0) then
      qtmx = 0._r8
      kmx = 0
      qtmn = 0._r8
      kmn = 0
      do k = 1,pver
         !        write (iulog,*) 'aaa ', k, ptend%q(icol,k,1), ptend%q(icol,k,2)
         if (ptend%q(icol,k,1).ge.qtmx) then
            kmx = k
            qtmx = ptend%q(icol,k,1)
         endif
         if (ptend%q(icol,k,1).le.qtmn) then
            kmn = k
            qtmn = ptend%q(icol,k,1)
         endif
      end do
      k = kmx
66    format ('tphysbc, aft shallow:', 4i5, 6f9.4)
      call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
      write (iulog,66)  nstep, icol, &
           kmx, kmn, &
           prec_cmf(icol)*8.64e7_r8, rliq2(icol)*8.64e7_r8,  &
           qtmx*8.64e7_r8, qtmn*8.64e7_r8, &
           (state%q(icol,k,1)+ptend%q(icol,k,1)*ztodt)*1.e3_r8, &
           (state%q(icol,k,ixcldliq)+ptend%q(icol,k,ixcldliq)*ztodt)*1.e3_r8

   endif

end subroutine phys_debug_shallow1

!================================================================================

subroutine phys_debug_strat1(state, ptend, nstep, prec_str, rliq, ztodt, kmx)

   use constituents, only: cnst_get_ind

   type(physics_state), intent(in) :: state             ! Physics state variables
   type(physics_ptend), intent(in) :: ptend             ! Physics process tendencies
   integer,             intent(in) :: nstep
   real(r8),            intent(in) :: prec_str(pcols)   ! sfc flux of precip from stratiform (m/s)
   real(r8),            intent(in) :: rliq(pcols)       ! vertical integral of liquid not yet in q(ixcldliq)
   real(r8),            intent(in) :: ztodt             ! physics time step
   integer,             intent(in) :: kmx

   integer  :: icol, k, ixcldliq
   !-----------------------------------------------------------------------------

   icol = phys_debug_col(state%lchnk)
   if (icol > 0) then

      k = kmx
67    format ('tphysbc, aft strat:', i5, 6f9.4)
      call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
      write (iulog,67)  nstep, prec_str(icol)*8.64e7_r8, rliq(icol)*8.64e7_r8,  &
           (ptend%q(icol,k,1)*ztodt)*1.e3_r8, &
           (ptend%q(icol,k,ixcldliq)*ztodt)*1.e3_r8, &
           (state%q(icol,k,1)+ptend%q(icol,k,1)*ztodt)*1.e3_r8, &
           (state%q(icol,k,ixcldliq)+ptend%q(icol,k,ixcldliq)*ztodt)*1.e3_r8

   endif

end subroutine phys_debug_strat1

!================================================================================

subroutine phys_debug_srf1(lchnk, cam_in)

   integer,             intent(in) :: lchnk             ! local chunk index
   type(cam_in_t),      intent(in) :: cam_in            ! CAM's import state

   integer :: icol
   !-----------------------------------------------------------------------------

   icol = phys_debug_col(lchnk)
   if (icol > 0) then

      write(iulog,*) 'bot tphysbc: cam_in%tref', cam_in%tref(icol)

   endif

end subroutine phys_debug_srf1

!================================================================================

subroutine phys_debug_hbdiff1(lchnk, pblh, zl, zh)

   integer,  intent(in) :: lchnk             ! local chunk index
   real(r8), intent(in) :: pblh(pcols)       ! boundary-layer height [m]
   real(r8), intent(in) :: zl(pcols)         ! zmzp / Obukhov length
   real(r8), intent(in) :: zh(pcols)         ! zmzp / pblh

   integer :: icol
   !-----------------------------------------------------------------------------

   icol = phys_debug_col(lchnk)
   if (icol > 0) then

      write(iulog,*) ' austach_pbl, pblh, zl, zh: ', pblh(icol), zl(icol), zh(icol)

   endif

end subroutine phys_debug_hbdiff1

!================================================================================

subroutine phys_debug_flux1(lchnk, srfflx, lhflx, shflx, taux, tauy, qflx, &
                            lhflx_res, shflx_res, taux_res, tauy_res, qflx_res)

   integer,            intent(in) :: lchnk      ! local chunk index
   type(cam_in_t),     intent(in) :: srfflx     ! cam import state
   real(r8),           intent(in) :: lhflx(:)   ! latent heat flux
   real(r8),           intent(in) :: shflx(:)   ! sensible heat flux
   real(r8),           intent(in) :: taux(:)    ! x momentum flux
   real(r8),           intent(in) :: tauy(:)    ! y momentum flux
   real(r8),           intent(in) :: qflx(:)    ! water vapor heat flux
   real(r8),           intent(in) :: lhflx_res(:)   ! latent heat flux
   real(r8),           intent(in) :: shflx_res(:)   ! sensible heat flux
   real(r8),           intent(in) :: taux_res(:)    ! x momentum flux
   real(r8),           intent(in) :: tauy_res(:)    ! y momentum flux
   real(r8),           intent(in) :: qflx_res(:)    ! water vapor heat flux

   integer :: icol
   !-----------------------------------------------------------------------------

   icol = phys_debug_col(lchnk)
   if (icol > 0) then

      write(iulog,*) ' b flux_tweak called, lhflx, oldlhflx ', &
          srfflx%lhf(icol), lhflx(icol)
      write(iulog,*) ' sfmodel fluxes lhf, shf, taux, tauy, q: ', &
          srfflx%lhf(icol), &
          srfflx%shf(icol), &
          srfflx%wsx(icol), &
          srfflx%wsy(icol), &
          srfflx%cflx(icol,1)
      write(iulog,*) ' last fluxes used lhf, shf, taux, tauy, q: ', &
          lhflx(icol), &
          shflx(icol), &
          taux(icol), &
          tauy(icol), &
          qflx(icol)
      write(iulog,*) ' current residuals lhf, shf, taux, tauy, q: ', &
          lhflx_res(icol), &
          shflx_res(icol), &
          taux_res(icol), &
          tauy_res(icol), &
          qflx_res(icol)

   endif

end subroutine phys_debug_flux1

!================================================================================

subroutine phys_debug_flux2(lchnk, srfflx, lhflx, &
                            lhflx_res, shflx_res, taux_res, tauy_res, qflx_res)

   integer,            intent(in) :: lchnk      ! local chunk index
   type(cam_in_t),     intent(in) :: srfflx     ! cam import state
   real(r8),           intent(in) :: lhflx(:)   ! latent heat flux
   real(r8),           intent(in) :: lhflx_res(:)   ! latent heat flux
   real(r8),           intent(in) :: shflx_res(:)   ! sensible heat flux
   real(r8),           intent(in) :: taux_res(:)    ! x momentum flux
   real(r8),           intent(in) :: tauy_res(:)    ! y momentum flux
   real(r8),           intent(in) :: qflx_res(:)    ! water vapor heat flux

   integer :: icol
   !-----------------------------------------------------------------------------

   icol = phys_debug_col(lchnk)
   if (icol > 0) then

      write(iulog,*) ' a flux_tweak called, lhflx, oldlhflx ', &
          srfflx%lhf(icol), lhflx(icol)
      write (iulog,66) ' residual fractions lhf, shf, taux, tauy, q ', &
          lhflx_res(icol)/srfflx%lhf(icol), &
          shflx_res(icol)/srfflx%shf(icol), &
          taux_res(icol)/srfflx%wsx(icol), &
          tauy_res(icol)/srfflx%wsy(icol), &
          qflx_res(icol)/srfflx%cflx(icol,1)
      write(iulog,66) ' residual lhf, shf, taux, tauy, q: ', &
          lhflx_res(icol), &
          shflx_res(icol), &
          taux_res(icol), &
          tauy_res(icol), &
          qflx_res(icol)
      write(iulog,*) ' used fluxes lhf, shf, taux, tauy, q: ', &
          srfflx%lhf(icol), &
          srfflx%shf(icol), &
          srfflx%wsx(icol), &
          srfflx%wsy(icol), &
          srfflx%cflx(icol,1)
66 format (a, 1p, 5e15.5)

   endif

end subroutine phys_debug_flux2

end module phys_debug
